home *** CD-ROM | disk | FTP | other *** search
/ Super PC 34 / Super PC 34 (Shareware).iso / spc / UTIL / DJGPP2 / V2 / DJTST200.ZIP / tests / libc / ansi / math / cephes / mtstold.c < prev   
Encoding:
C/C++ Source or Header  |  1993-08-08  |  10.4 KB  |  431 lines
  1. /*   mtst.c
  2.  Consistency tests for math functions.
  3.  Inverse function tests and Wronksians for random arguments.
  4.  With NTRIALS=10000, the following are typical results for
  5.  IEEE double precision arithmetic:
  6.  
  7. Consistency test of math functions.
  8. Max and rms relative errors for 10000 random arguments.
  9. x =   cbrt(   cube(x) ):  max = 1.54E-016   rms = 6.50E-018
  10. x =   atan(    tan(x) ):  max = 2.83E-016   rms = 7.68E-017
  11. x =    sin(   asin(x) ):  max = 4.44E-016   rms = 7.01E-017
  12. x =   sqrt( square(x) ):  max = 1.57E-016   rms = 2.70E-017
  13. x =    log(    exp(x) ):  max = 1.88E-014   rms = 2.02E-016
  14. x =   tanh(  atanh(x) ):  max = 4.00E-016   rms = 8.29E-017
  15. x =  asinh(   sinh(x) ):  max = 2.04E-016   rms = 1.56E-017
  16. x =  acosh(   cosh(x) ):  max = 2.67E-014   rms = 2.95E-016
  17. x = pow( pow(x,a),1/a ):  max = 3.38E-012   rms = 3.39E-014
  18. Legendre  ellpk,  ellpe:  max = 1.55E-011   rms = 1.74E-013
  19. x =  ndtri(   ndtr(x) ):  max = 4.09E-014   rms = 6.40E-016
  20. Absolute error criterion (but relative if >1):
  21. lgam(x) = log(gamma(x)):  max = 4.49E-016   rms = 8.80E-017
  22. x =    cos(   acos(x) ):  max = 4.44E-016   rms = 1.00E-016
  23. Absolute error and only 2000 trials:
  24. Wronksian of   Yn,   Jn:  max = 7.12E-011   rms = 1.59E-012
  25. Wronksian of   Kn,   Iv:  max = 7.07E-012   rms = 5.25E-013
  26.  
  27. */
  28.  
  29. /*
  30. Cephes Math Library Release 2.1:  December, 1988
  31. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  32. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  33. f*/
  34.  
  35. #include <stdlib.h>
  36. #include <stdio.h>
  37. #include <math.h>
  38. #include <string.h>
  39. #include <signal.h>
  40. #include <mathext.h>
  41. #include "mconf.h"
  42.  
  43. #define NTRIALS 10000
  44. #define WTRIALS (NTRIALS/5)
  45. #define STRTST 0
  46.  
  47. double          ndtr(), ndtri(), ellpe(), ellpk(), gamma();
  48. double          fabs(), sqrt();
  49. double          cbrt(), exp(), log(), tan(), atan();
  50. double          sin(), asin(), cos(), acos(), pow();
  51. double          tanh(), atanh(), sinh(), asinh(), cosh(), acosh();
  52. double          jn(), yn(), iv(), kn();
  53.  
  54. /* Provide inverses for square root and cube root: */
  55. double 
  56. square(x)
  57.     double          x;
  58. {
  59.     return (x * x);
  60. }
  61.  
  62. double
  63. cube(x)
  64.     double          x;
  65. {
  66.     return (x * x * x);
  67. }
  68.  
  69. /* lookup table for each function */
  70. struct fundef
  71. {
  72.     char           *nam1;    /* the function */
  73.     double          (*name) ();
  74.     char           *nam2;    /* its inverse  */
  75.     double          (*inv) ();
  76.     int             nargs;    /* number of function arguments */
  77.     int             tstyp;    /* type code of the function */
  78.     long            ctrl;    /* relative error flag */
  79.     double          arg1w;    /* width of domain for 1st arg */
  80.     double          arg1l;    /* lower bound domain 1st arg */
  81.     long            arg1f;    /* flags, e.g. integer arg */
  82.     double          arg2w;    /* same info for args 2, 3, 4 */
  83.     double          arg2l;
  84.     long            arg2f;
  85. /*
  86.     double arg3w;
  87.     double arg3l;
  88.     long arg3f;
  89.     double arg4w;
  90.     double arg4l;
  91.     long arg4f;
  92. */
  93. };
  94.  
  95.  
  96. /* fundef.ctrl bits: */
  97. #define RELERR 1
  98. #define EXPSCAL 4
  99.  
  100. /* fundef.tstyp  test types: */
  101. #define POWER 1
  102. #define ELLIP 2
  103. #define GAMMA 3
  104. #define WRONK1 4
  105. #define WRONK2 5
  106. #define WRONK3 6
  107.  
  108. /* fundef.argNf  argument flag bits: */
  109. #define INT 2
  110.  
  111. extern double   MINLOG;
  112. extern double   MAXLOG;
  113. extern double   CPI;
  114. extern double   PIO2;
  115.  
  116. /*
  117. define MINLOG -170.0
  118. define MAXLOG +170.0
  119. define PI 3.14159265358979323846
  120. define PIO2 1.570796326794896619
  121. */
  122.  
  123. #define NTESTS 12
  124. struct fundef   defs[NTESTS] =
  125. {
  126.  {"  cube",   cube,   "  cbrt",   cbrt, 1, 0, 1, 2000.0, -1000.0,   0,
  127. 0.0, 0.0, 0},
  128.     {"   tan", tan, "  atan", atan, 1, 0, 1, 0.0, 0.0, 0,
  129.      0.0, 0.0, 0},
  130.     {"  asin", asin, "   sin", sin, 1, 0, 1, 2.0, -1.0, 0,
  131.      0.0, 0.0, 0},
  132.     {"square", square, "  sqrt", sqrt, 1, 0, 1, 170.0, -85.0, EXPSCAL,
  133.      0.0, 0.0, 0},
  134.     {"   exp", exp, "   log", log, 1, 0, 1, 340.0, -170.0, 0,
  135.      0.0, 0.0, 0},
  136.     {" atanh", atanh, "  tanh", tanh, 1, 0, 1, 2.0, -1.0, 0,
  137.      0.0, 0.0, 0},
  138.     {"  sinh", sinh, " asinh", asinh, 1, 0, 1, 340.0, 0.0, 0,
  139.      0.0, 0.0, 0},
  140.     {"  cosh", cosh, " acosh", acosh, 1, 0, 1, 340.0, 0.0, 0,
  141.      0.0, 0.0, 0},
  142.     {"pow", pow, "pow", pow, 2, POWER, 1, 25.0, 0.0, 0,
  143.      50.0, -25.0, 0},
  144. /* {" ellpe",  ellpe,   " ellpk",  ellpk, 1, ELLIP, 1, 1.0,    0.0,  0,
  145. 0.0, 0.0, 0},
  146. {"  ndtr",   ndtr,   " ndtri",  ndtri, 1, 0, 1,   10.0,   -10.0,   0,
  147. 0.0, 0.0, 0}, */
  148. /* Absolute error criterion starts with tstyp = GAMMA: */
  149.     {"gamma", gamma, "lgamma", lgamma, 1, GAMMA, 0, 11.0, 1.0, 0,
  150.      0.0, 0.0, 0},
  151.     {"  acos", acos, "   cos", cos, 1, 0, 0, 2.0, -1.0, 0,
  152.      0.0, 0.0, 0},
  153. /* Smaller number of trials starts with tstyp = WRONK1: */
  154.     {"  Jnt", jn, "  Ynt", yn, 2, WRONK1, 0, 30.0, 0.1, 0,
  155.      /* 40.0, -20.0,*/ 20.0, 0.0, INT},
  156. /* {"  Iv",     iv,   "  Kn",     kn, 2, WRONK2, 0,  30.0, 0.1,  0,
  157. 20.0, 0.0, INT}, */
  158. };
  159.  
  160. static char    *headrs[] =
  161. {
  162.     "x = %s( %s(x) ): ",
  163.     "x = %s( %s(x,a),1/a ): ",    /* power */
  164.     "Legendre %s, %s: ",    /* ellip */
  165.     "%s(x) = log(%s(x)): ",    /* gamma */
  166.     "Wronksian of %s, %s: ",
  167.     "Wronksian of %s, %s: ",
  168.     "Wronksian of %s, %s: "
  169. };
  170.  
  171. static double   cy1 = 0.0;
  172. static double   y2 = 0.0;
  173. static double   y3 = 0.0;
  174. static double   y4 = 0.0;
  175. static double   a = 0.0;
  176. static double   b = 0.0;
  177. static double   c = 0.0;
  178. static double   x = 0.0;
  179. static double   y = 0.0;
  180. static double   z = 0.0;
  181. static double   e = 0.0;
  182. static double   max = 0.0;
  183. static double   rmsa = 0.0;
  184. static double   rms = 0.0;
  185. static double   ave = 0.0;
  186.  
  187. void mysig(int sig);
  188. int
  189. main()
  190. {
  191.     double          (*fun) ();
  192.     double          (*ifun) ();
  193.     struct fundef  *d;
  194.     int             i, j, k, itst;
  195.     int             m, ntr;
  196.     _set_cw87(0x137f);
  197.     ntr = NTRIALS;
  198.     printf("Consistency test of math functions.\n");
  199.     printf("Max and rms relative errors for %d random arguments.\n",
  200.        ntr);
  201.  
  202. /* Initialize machine dependent parameters: */
  203. #ifdef __TURBOC__
  204. #define FACT 0.99
  205. #else
  206. #define FACT 1.0000
  207. #endif
  208.     defs[1].arg1w = CPI;
  209.     defs[1].arg1l = -CPI/2.0;
  210.     defs[3].arg1w = MAXLOG*FACT;
  211.     defs[3].arg1l = -MAXLOG/2.0*FACT;
  212.     defs[4].arg1w = 2.0*MAXLOG*FACT;
  213.     defs[4].arg1l = -MAXLOG*FACT;
  214.     defs[6].arg1w = 2.0*MAXLOG*FACT;
  215.     defs[6].arg1l = -MAXLOG*FACT;
  216.     defs[7].arg1w = MAXLOG*FACT;
  217.     defs[7].arg1l = 0.0;
  218.  
  219.  
  220. /* Outer loop, on the test number: */
  221.  
  222.     for (itst = STRTST ; itst < NTESTS; itst++)
  223.     {
  224.     d = &defs[itst];
  225.     m = 0;
  226.     max = 0.0;
  227.     rmsa = 0.0;
  228.     ave = 0.0;
  229.     fun = d->name;
  230.     ifun = d->inv;
  231.  
  232. /* Absolute error criterion starts with gamma function
  233.  * (put all such at end of table)
  234.  */
  235.     if (d->tstyp == GAMMA)
  236.         printf("Absolute error criterion (but relative if >1):\n");
  237.  
  238. /* Smaller number of trials for Wronksians
  239.  * (put them at end of list)
  240.  */
  241.     if (d->tstyp == WRONK1)
  242.     {
  243.         ntr = WTRIALS;
  244.         printf("Absolute error and only %d trials:\n", ntr);
  245.     }
  246.  
  247.     printf(headrs[d->tstyp], d->nam2, d->nam1);
  248.     printf("random arguments in %g %g\n",
  249.         d->arg1l, d->arg1l + d->arg1w);
  250.     for (i = 0; i < ntr; i++)
  251.     {
  252.         m++;
  253.  
  254. /* make random number(s) in desired range(s) */
  255.         switch (d->nargs)
  256.         {
  257.  
  258.         default:
  259.             goto illegn;
  260.  
  261.         case 2:
  262.             drand(&a);
  263.             a = d->arg2w * (a - 1.0) + d->arg2l;
  264.             if (d->arg2f & EXPSCAL)
  265.             {
  266.             a = exp(a);
  267.             drand(&y2);
  268.             a -= 1.0e-13 * a * y2;
  269.             }
  270.             if (d->arg2f & INT)
  271.             {
  272.             k = a + 0.25;
  273.             a = k;
  274.             }
  275.  
  276.         case 1:
  277.             drand(&x);
  278.             x = d->arg1w * (x - 1.0) + d->arg1l;
  279.             if (d->arg1f & EXPSCAL)
  280.             {
  281.             x = exp(x);
  282.             drand(&a);
  283.             x += 1.0e-13 * x * a;
  284.             }
  285.         }
  286.  
  287.  
  288. /* compute function under test */
  289.         switch (d->nargs)
  290.         {
  291.         case 1:
  292.             switch (d->tstyp)
  293.             {
  294.             case ELLIP:
  295.                 cy1 = (*(fun)) (x);
  296.                 y2 = (*(fun)) (1.0 - x);
  297.                 y3 = (*(ifun)) (x);
  298.                 y4 = (*(ifun)) (1.0 - x);
  299.                 break;
  300.  
  301.             case GAMMA:
  302.                 y = lgamma(x);
  303.                 x = log(gamma(x));
  304.                 break;
  305.  
  306.             default:
  307.                 z = (*(fun)) (x);
  308.                 y = (*(ifun)) (z);
  309.             }
  310.             break;
  311.  
  312.         case 2:
  313.             if (d->arg2f & INT)
  314.             {
  315.             switch (d->tstyp)
  316.             {
  317.                 case WRONK1:
  318.                 cy1 = (*fun) (k, x);    /* jn */
  319.                 y2 = (*fun) (k + 1, x);
  320.                 y3 = (*ifun) (k, x);    /* yn */
  321.                 y4 = (*ifun) (k + 1, x);
  322.                 break;
  323.  
  324.                 case WRONK2:
  325.                 cy1 = (*fun) (a, x);    /* iv */
  326.                 y2 = (*fun) (a + 1.0, x);
  327.                 y3 = (*ifun) (k, x);    /* kn */
  328.                 y4 = (*ifun) (k + 1, x);
  329.                 break;
  330.  
  331.                 default:
  332.                 z = (*fun) (k, x);
  333.                 y = (*ifun) (k, z);
  334.             }
  335.             }
  336.             else
  337.             {
  338.             if (d->tstyp == POWER)
  339.             {
  340.                 z = (*fun) (x, a);
  341.                 y = (*ifun) (z, 1.0 / a);
  342.             }
  343.             else
  344.             {
  345.                 z = (*fun) (a, x);
  346.                 y = (*ifun) (a, z);
  347.             }
  348.             }
  349.             break;
  350.  
  351.  
  352.         default:
  353.           illegn:
  354.             printf("Illegal nargs= %d", d->nargs);
  355.             exit(1);
  356.         }
  357.  
  358.         switch (d->tstyp)
  359.         {
  360.         case WRONK1:
  361.             e = (y2 * y3 - cy1 * y4) - 2.0 / (CPI * x);    /* Jn, Yn */
  362.             break;
  363.  
  364.         case WRONK2:
  365.             e = (y2 * y3 + cy1 * y4) - 1.0 / x;    /* In, Kn */
  366.             break;
  367.  
  368.         case ELLIP:
  369.             e = (cy1 - y3) * y4 + y3 * y2 - PIO2;
  370.             break;
  371.  
  372.         default:
  373.             e = y - x;
  374.             break;
  375.         }
  376.  
  377.         if (d->ctrl & RELERR)
  378.         e /= x;
  379.         else
  380.         {
  381.         if (fabs(x) > 1.0)
  382.             e /= x;
  383.         }
  384.  
  385.         ave += e;
  386. /* absolute value of error */
  387.         if (e < 0)
  388.         e = -e;
  389.  
  390. /* peak detect the error */
  391.         if (e > max)
  392.         {
  393.         max = e;
  394.  
  395.         if (e > 1.0e-10)
  396.         {
  397.             /* printf("x %.6E z %.6E y %.6E max %.4E\n",
  398.                x, z, y, max); */
  399.             printf("x %.18E z %.18E y %.18E max %.4E\n",
  400.                x, z, y, max);
  401.             if (d->tstyp >= WRONK1)
  402.             {
  403.             printf("cy1 %.4E y2 %.4E y3 %.4E y4 %.4E k %d x %.4E\n",
  404.                    cy1, y2, y3, y4, k, x);
  405.             }
  406.         }
  407.  
  408. /*
  409.     printf("%.8E %.8E %.4E %6ld \n", x, y, max, n);
  410.     printf("%d %.8E %.8E %.4E %6ld \n", k, x, y, max, n);
  411.     printf("%.6E %.6E %.6E %.4E %6ld \n", a, x, y, max, n);
  412.     printf("%.6E %.6E %.6E %.6E %.4E %6ld \n", a, b, x, y, max, n);
  413.     printf("%.4E %.4E %.4E %.4E %.4E %.4E %6ld \n",
  414.         a, b, c, x, y, max, n);
  415. */
  416.         }
  417.  
  418. /* accumulate rms error    */
  419.         e *= 1.0e16;    /* adjust range */
  420.         rmsa += e * e;    /* accumulate the square of the error */
  421.       endlup:
  422.         ;
  423.     }
  424.  
  425. /* report after NTRIALS trials */
  426.     rms = 1.0e-16 * sqrt(rmsa / m);
  427.     printf(" max = %.2E   rms = %.2E\n", max, rms);
  428.     }                /* loop on itst */
  429.     return(0);
  430. }
  431.